Thermal performance: CFUs and growth curves time-to-detection analyses
Hermina and Anjaney
Introduction
This notebook is for analyzing the temperature growth curves data
produced by Anjaney and Hermina, as well as the survival CFU data
produced by Hermina. It summarizes the findings for publication. This
notebook was created by copying the relevant sections from
Anjaney_analysis--hermina_redo.Rmd and
plate_reader_calibration.Rmd.
The goal is to describe the thermal niche of different Pseudomonas strains present in the Terrestrial Ecology lab of Maddy Thakur (University of Bern). The main question the experiments are trying to answer: Is there a relationship between a strain’s growth (fast VS slow grower) and its thermal niche?
See the associated manuscript for more details. This notebook covers all analyses performed for Experiment I in the manuscript.
The data
16 soil Pseudomonas strains were used in total (BSC001 - BSC010, BSC015, BSC019, CK101 - CK104), inoculated either from early exponential phase or from stationary phase.
There are 3 data sets which are integrated in this analysis:
Colony forming unit (CFU) data for all 16 strains growing at 28C, 35C, and 40C.
Growth curves inoculated in a 10-fold serial dilution at 6 different starting concentrations (\(N_0\)) and 4 different temperatures (25C, 30C, 35C, and 40C). The CFU data is integrated into this analysis because it is used to infer the OD-converted inoculum size for each strain x temperature combination.
Dilution series performed on both of the microspectrophotometers in our lab. This is used to integrate the CFU data with the OD data (see above).
All data files and protocols are associated with the above analyses are available in the folder and sub-folders where this file is found.
Before we get started, let’s load all the R packages and list the versions for reproducibility of the entire environment:
# load environment
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(readxl) # for importing data directly from Excel sheet
library(broom) # for easy fitting of several lm's from the same df
library(Hmisc) # for calculating non-parametric bootstrap confidence intervals##
## Attaching package: 'Hmisc'
##
## The following objects are masked from 'package:dplyr':
##
## src, summarize
##
## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(AICcmodavg) # AIC comparison for survival CFU data
library("CarletonStats") # for bootstrapping the slope of linear regression
library(DescTools) # for Kendall's Coefficient of Concordance (W)##
## Attaching package: 'DescTools'
##
## The following objects are masked from 'package:Hmisc':
##
## %nin%, Label, Mean, Quantile# this is a set of custom functions written by AHG:
source("H1_extract_data_fun_v3-1.R") # for getting the data into R from the default H1 txt export format
# print the complete info about packages and versions currently loaded in the environment:
sessionInfo()## R version 4.4.2 (2024-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_United Kingdom.utf8
## [2] LC_CTYPE=English_United Kingdom.utf8
## [3] LC_MONETARY=English_United Kingdom.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United Kingdom.utf8
##
## time zone: Europe/Paris
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] DescTools_0.99.58 CarletonStats_2.2 AICcmodavg_2.3-3 Hmisc_5.2-1
## [5] broom_1.0.7 readxl_1.4.3 lubridate_1.9.3 forcats_1.0.0
## [9] stringr_1.5.1 dplyr_1.1.4 purrr_1.0.2 readr_2.1.5
## [13] tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] Exact_3.3 tidyselect_1.2.1 rootSolve_1.8.2.4 farver_2.1.2
## [5] fastmap_1.2.0 digest_0.6.37 rpart_4.1.23 timechange_0.3.0
## [9] lifecycle_1.0.4 cluster_2.1.6 survival_3.7-0 lmom_3.2
## [13] magrittr_2.0.3 compiler_4.4.2 rlang_1.1.4 sass_0.4.9
## [17] tools_4.4.2 yaml_2.3.10 data.table_1.16.2 knitr_1.49
## [21] htmlwidgets_1.6.4 expm_1.0-0 withr_3.0.2 foreign_0.8-87
## [25] nnet_7.3-19 grid_4.4.2 stats4_4.4.2 unmarked_1.4.3
## [29] xtable_1.8-4 e1071_1.7-16 colorspace_2.1-1 scales_1.3.0
## [33] MASS_7.3-61 cli_3.6.3 mvtnorm_1.3-2 rmarkdown_2.29
## [37] generics_0.1.3 rstudioapi_0.17.1 httr_1.4.7 tzdb_0.4.0
## [41] gld_2.6.6 cachem_1.1.0 proxy_0.4-27 splines_4.4.2
## [45] parallel_4.4.2 cellranger_1.1.0 base64enc_0.1-3 vctrs_0.6.5
## [49] boot_1.3-31 Matrix_1.7-1 jsonlite_1.8.9 VGAM_1.1-12
## [53] hms_1.1.3 patchwork_1.3.0 Formula_1.2-5 htmlTable_2.4.3
## [57] jquerylib_0.1.4 glue_1.8.0 stringi_1.8.4 gtable_0.3.6
## [61] munsell_0.5.1 pillar_1.10.1 htmltools_0.5.8.1 R6_2.5.1
## [65] evaluate_1.0.1 lattice_0.22-6 haven_2.5.4 backports_1.5.0
## [69] bslib_0.8.0 class_7.3-22 Rcpp_1.0.13-1 gridExtra_2.3
## [73] nlme_3.1-166 checkmate_2.3.2 xfun_0.49 pkgconfig_2.0.3# set theme for all plots
fave_theme <- theme_light() + # see other options at https://ggplot2.tidyverse.org/reference/ggtheme.html
theme(text = element_text(size=15), # larger text size for titles & axes
panel.grid.major = element_blank(), # remove major gridlines
panel.grid.minor = element_blank()) # remove minor gridlines
theme_set(fave_theme)
# define a palette for plotting all 6 species
ALLsp_palette = palette.colors(8, palette = "R4")[2:7]Plate reader calibration for CFU-to-OD conversion
Analysis of data generated by following the protocol
OD Calibration to CFUs.docx.
Fairly exhaustive calibrations were done for the microplate spectrophotometers (hereafter referred to generally as “plate readers” or specifically by the name of each instrument, BioTek Synergy H1 or Biotek Epoch2). Different Pseudomonas species were grown in different culture conditions (96-well microplate, tube, or Erlenmeyer), temperatures (25C or 40C), and culture phases (exponential or stationary). Then, dilution series were made. These were measured by spectrophotometer at 600nm (OD\(_{600}\)) and by CFU to get the calibration. There were a lot of inconsistencies; these are to be expected since OD measures the amount of scattered light, which depends in part on cell density but also on cell size, clumping, etc.
Anyway, I tried my best to estimate a calibration curve.
Load & parse data
#####################################
# load the microplate OD data
#####################################
synergy.df <- read_excel("./data/calibration_spectrophotometers--22Feb24.xlsx", sheet="microplate H1") %>%
pivot_longer(cols=`25_rep1`:`40_rep6`,
names_to=c("Temp", "Rep"),
names_pattern="(.*)_rep(.)",
values_to="OD") # reorganize the shape of the dataframe
epoch.df <- read_excel("./data/calibration_spectrophotometers--22Feb24.xlsx", sheet="microplate Epoch") %>%
pivot_longer(cols=`25_rep1`:`40_rep6`,
names_to=c("Temp", "Rep"),
names_pattern="(.*)_rep(.)",
values_to="OD") # reorganize the shape of the dataframe
# load the microplate CFU data
plateCFU.df <- read_excel("./data/calibration_spectrophotometers--22Feb24.xlsx", sheet="microplate_CFU")
#####################################
# baseline subtract using the blanks as well as the bottom 3 lowest dilutions for each sample
#####################################
## check if there's any systematic trend in the baseline value
#baseline_trend <- rbind(synergy.df %>% add_column(machine="Synergy"), # combine the data from the 2 plate readers
# epoch.df %>% add_column(machine="Epoch")) %>%
# filter(Dilution < 5e-06) %>% # keep blanks & 3 lowest dilutions
# # get median value for each temperature & replicate
# select(-Dilution, -Sample) %>% group_by(Temp, Rep, machine) %>%
# summarise(medianBlankOD = median(OD, na.rm=TRUE),
# ICR_lo = quantile(OD, probs=0.25),
# ICR_hi = quantile(OD, probs=0.75))
#ggplot(baseline_trend, aes(x=Rep, y=medianBlankOD, colour=Temp)) +
# facet_wrap(vars(machine), scales="free") +
# geom_errorbar(aes(ymin=ICR_lo, ymax=ICR_hi), width=0.2) +
# geom_point()
## there is NO systematic trend in the baseline value
#rm(baseline_trend)
# in this case it seems there's a difference in the baselines for different incubation temperatures, but only for the Epoch2...
# baseline everything in the same way as before
synergy.df <- inner_join(synergy.df %>% filter(Dilution < 5e-06) %>%
select(-Dilution, -Sample) %>% # get the blanks & calculate mean value for each Temp & Rep
group_by(Temp, Rep) %>% summarise(medianBlankOD = median(OD, na.rm=TRUE)),
synergy.df %>% filter(Dilution > 5e-06) # add the medianBlankOD back to the data without the samples used for blanking
) %>% mutate(baselinedOD = OD-medianBlankOD)## `summarise()` has grouped output by 'Temp'. You can override using the
## `.groups` argument.
## Joining with `by = join_by(Temp, Rep)`epoch.df <- inner_join(epoch.df %>% filter(Dilution < 5e-06) %>%
select(-Dilution, -Sample) %>% # get the blanks & calculate mean value for each Temp & Rep
group_by(Temp, Rep) %>%
summarise(medianBlankOD = median(OD, na.rm=TRUE)),
epoch.df %>% filter(Dilution > 5e-06) # add the medianBlankOD back to the data without the samples used for blanking
) %>% mutate(baselinedOD = OD-medianBlankOD)## `summarise()` has grouped output by 'Temp'. You can override using the
## `.groups` argument.
## Joining with `by = join_by(Temp, Rep)`#####################################
# estimate the mean CFU of each sample
#####################################
plateCFU.df <- plateCFU.df %>% mutate(CFU_per_mL = `CFU per 100uL`/Dilution*10) %>%
group_by(Sample, Phase, Incubation) %>%
summarise(meanCFU = mean(CFU_per_mL),
sdCFU = sd(CFU_per_mL))## `summarise()` has grouped output by 'Sample', 'Phase'. You can override using
## the `.groups` argument.# combine the CFU data with the OD data & infer the CFU for each dilution used in the OD data
synergy.df <- left_join(synergy.df, plateCFU.df) %>%
mutate(inferred_meanDensity = meanCFU*Dilution,
inferred_sdDensity = sdCFU*Dilution) %>%
select(-meanCFU, -sdCFU)## Joining with `by = join_by(Sample, Incubation)`epoch.df <- left_join(epoch.df, plateCFU.df) %>%
mutate(inferred_meanDensity = meanCFU*Dilution,
inferred_sdDensity = sdCFU*Dilution) %>%
select(-meanCFU, -sdCFU)## Joining with `by = join_by(Sample, Incubation)`# remove some of the data points because it seems we are hitting up against the limit of detection
synergy.df <- synergy.df %>% filter(inferred_meanDensity > 5e+04,
baselinedOD > 5e-03,
Temp==25) # drop the 40C because it looks even worse than 25C
# remove some of the data points because it seems we are hitting up against the limit of detection
epoch.df <- epoch.df %>% filter(inferred_meanDensity > 1e+04,
baselinedOD > 5e-03,
Temp==25) # again remove OD acquisition temp of 40C
# create columns for log10 values
synergy.df <- synergy.df %>%
mutate(log10_OD = log10(baselinedOD),
log10_density = log10(inferred_meanDensity)) %>%
filter(log10_OD > -Inf)
epoch.df <- epoch.df %>%
mutate(log10_OD = log10(baselinedOD),
log10_density = log10(inferred_meanDensity)) %>%
filter(log10_OD > -Inf)
# store the raw data
all_calib.data <- rbind(synergy.df %>% mutate(Detector="synergy"),
epoch.df %>% mutate(Detector="epoch")) %>% ungroup() %>%
select(-Rep, -medianBlankOD, -Dilution, -OD, -Temp) %>%
separate_wider_delim(Incubation, "C ", names=c("Incubation_temp", "Flask")) %>%
mutate(OD_temp=25)
all_calib.data$Incubation_temp <- as.numeric(all_calib.data$Incubation_temp)
rm(epoch.df, synergy.df, plateCFU.df)Load additional data…
# load the microplate OD data
synergy.df <- read_excel("./data/calibration_spectrophotometers--26Feb24.xlsx", sheet="microplate H1") %>%
pivot_longer(cols =`rep1`:`rep6`,
names_to = "Rep",
names_pattern="rep(.)",
values_to="OD") %>% # reorganize the shape of the dataframe
mutate(OD_temp = 25)
epoch.df <- read_excel("./data/calibration_spectrophotometers--26Feb24.xlsx", sheet="microplate Epoch") %>%
pivot_longer(cols =`rep1`:`rep6`,
names_to = "Rep",
names_pattern="rep(.)",
values_to="OD") %>% # reorganize the shape of the dataframe
mutate(OD_temp = 40)
# load the microplate CFU data
plateCFU.df <- read_excel("./data/calibration_spectrophotometers--26Feb24.xlsx", sheet="microplate_CFU")
# baseline subtract using the blanks as well as the bottom 1 lowest dilution for each sample
## check if there's any systematic trend in the baseline value
#baseline_trend <- rbind(synergy.df %>% add_column(machine="Synergy"), # combine the data from the 2 plate readers
# epoch.df %>% add_column(machine="Epoch")) %>%
# filter(Dilution < 1e-06) %>% # keep blanks & 3 lowest dilutions
# # get median value for each reader & replicate
# select(-Dilution, -Sample) %>% group_by(Rep, machine, Incubation) %>%
# summarise(medianBlankOD = median(OD, na.rm=TRUE),
# ICR_lo = quantile(OD, probs=0.25, na.rm=TRUE),
# ICR_hi = quantile(OD, probs=0.75, na.rm=TRUE))
#ggplot(baseline_trend, aes(x=Rep, y=medianBlankOD, colour=Incubation)) +
# facet_wrap(vars(machine), scales="free") +
# geom_errorbar(aes(ymin=ICR_lo, ymax=ICR_hi), width=0.2) +
# geom_point()
## yes, there is a trend. But only for the Epoch2 plate reader and 25C data...
#rm(baseline_trend)
# remove the first 5 replicates in the "25C plate" Epoch data
epoch.df <- rbind(epoch.df %>% filter(Incubation == "25C plate reader", Rep==6),
epoch.df %>% filter(Incubation == "40C plate reader", Rep<6))
# & remove the 6th replicate in the "40C plate" Epoch data (that one is all NA's anyway)
# now, after excluding problematic values, baseline similarly to previous data:
synergy.df <- inner_join(synergy.df %>% filter(Dilution < 1e-06) %>%
select(-Dilution, -Sample) %>% # get the blanks & calculate mean value for each Temp & Rep
group_by(Rep, Incubation) %>% summarise(medianBlankOD = median(OD, na.rm=TRUE)),
synergy.df %>% filter(Dilution > 1e-06) # add the medianBlankOD back to the data without the samples used for blanking
) %>% mutate(baselinedOD = OD-medianBlankOD)## `summarise()` has grouped output by 'Rep'. You can override using the `.groups`
## argument.
## Joining with `by = join_by(Rep, Incubation)`epoch.df <- inner_join(epoch.df %>% filter(Dilution < 1e-06) %>%
select(-Dilution, -Sample) %>% # get the blanks & calculate mean value for each Temp & Rep
group_by(Rep, Incubation) %>% summarise(medianBlankOD = median(OD, na.rm=TRUE)),
epoch.df %>% filter(Dilution > 1e-06) # add the medianBlankOD back to the data without the samples used for blanking
) %>% mutate(baselinedOD = OD-medianBlankOD)## `summarise()` has grouped output by 'Rep'. You can override using the `.groups`
## argument.
## Joining with `by = join_by(Rep, Incubation)`# estimate the mean CFU of each sample
plateCFU.df <- plateCFU.df %>% mutate(CFU_per_mL = `CFU per 100uL`/Dilution*10) %>%
group_by(Sample, Phase, Incubation) %>%
summarise(meanCFU = mean(CFU_per_mL),
sdCFU = sd(CFU_per_mL))## `summarise()` has grouped output by 'Sample', 'Phase'. You can override using
## the `.groups` argument.# combine the CFU data with the OD data & infer the CFU for each dilution used in the OD data
synergy.df <- left_join(synergy.df, plateCFU.df) %>%
mutate(inferred_meanDensity = meanCFU*Dilution,
inferred_sdDensity = sdCFU*Dilution) %>%
select(-meanCFU, -sdCFU)## Joining with `by = join_by(Incubation, Sample, Phase)`epoch.df <- left_join(epoch.df, plateCFU.df) %>%
mutate(inferred_meanDensity = meanCFU*Dilution,
inferred_sdDensity = sdCFU*Dilution) %>%
select(-meanCFU, -sdCFU)## Joining with `by = join_by(Incubation, Sample, Phase)`# remove some of the data points because it seems we are hitting up against the limit of detection
synergy.df <- rbind(synergy.df %>% filter(Sample=="BSC004", baselinedOD > 1e-02),
synergy.df %>% filter(Sample=="BSC008", Incubation=="25C plate reader", baselinedOD > 2e-02),
synergy.df %>% filter(Sample=="BSC008", Incubation=="40C plate reader", baselinedOD > 5e-02),
synergy.df %>% filter(Sample=="BSC015", Incubation=="25C plate reader", baselinedOD > 2e-02),
synergy.df %>% filter(Sample=="BSC015", Incubation=="40C plate reader", baselinedOD > 1e-02),
synergy.df %>% filter(Sample=="CK101", baselinedOD > 2e-02))
# remove the most concentrated samples because of violation of the Beer-Lambert law
synergy.df <- synergy.df %>% filter(baselinedOD < 1.4)
# remove some of the data points because it seems we are hitting up against the limit of detection
epoch.df <- rbind(epoch.df %>% filter(Sample %in% c("BSC004", "BSC008"), Incubation=="25C plate reader", baselinedOD > 5e-03),
epoch.df %>% filter(Sample %in% c("BSC015", "CK101"), Incubation=="25C plate reader", baselinedOD > 1e-02),
epoch.df %>% filter(Sample=="BSC015", Incubation=="40C plate reader", baselinedOD > 1e-02),
epoch.df %>% filter(Sample=="BSC008", Incubation=="40C plate reader", baselinedOD > 2e-02))
# remove the most concentrated samples because of violation of the Beer-Lambert law
epoch.df <- epoch.df %>% filter(baselinedOD < 1.1)
# get the log10 of the data
synergy.df <- synergy.df %>%
mutate(log10_OD = log10(baselinedOD),
log10_density = log10(inferred_meanDensity)) %>%
filter(log10_OD > -Inf,
baselinedOD > 1e-13)
epoch.df <- epoch.df %>%
mutate(log10_OD = log10(baselinedOD),
log10_density = log10(inferred_meanDensity)) %>%
filter(log10_OD > -Inf)
# store the raw data
temp.data <- rbind(synergy.df %>% mutate(Detector="synergy", OD_temp=25),
epoch.df %>% mutate(Detector="epoch", OD_temp=40)) %>% ungroup() %>%
select(-Rep, -medianBlankOD, -Dilution, -OD) %>%
separate_wider_delim(Incubation, "C ", names=c("Incubation_temp", "Flask"))
temp.data$Incubation_temp <- as.numeric(temp.data$Incubation_temp)
all_calib.data <- rbind(all_calib.data %>% mutate(Sample = "BSC015"),
temp.data)
rm(epoch.df, synergy.df, plateCFU.df, temp.data)Estimate CFU to OD conversion
Let’s consider only bacteria grown in the microplate because this is most relevant to our experiments. Then, we’re interested in averaging across the effects of temperature, species, and microplate reader. So let’s simply use data from different temperatures, species, and readers to run a linear model where CFU is the independent variable and OD is the dependent.
ggplot(all_calib.data %>% filter(Flask == "plate reader"),
aes(x=inferred_meanDensity, y=baselinedOD, colour=Detector)) +
geom_point() +
scale_y_log10() +
scale_x_log10() +
labs(title = "inocula grown in microplate")
## rather arbitrarily pick a subset of the data for creating the conversion
conversion.data <- all_calib.data %>% filter(Flask == "plate reader") %>%
filter(Sample != "BSC004") # exclude right side outlier
# exclude left side outlier
conversion.data <- conversion.data[-which(conversion.data$Sample=="BSC015" & conversion.data$Incubation_temp==40 & conversion.data$Phase=="Stationary"),]
# plot the final data that will be used for the conversion estimate:
ggplot(conversion.data, # here in terms of species
aes(x=inferred_meanDensity, y=baselinedOD, colour=Sample)) +
geom_point() +
scale_y_log10() +
scale_x_log10() +
labs(title = "data used for CFU to OD conversion")
# plot again in terms of the microplate reader
ggplot(conversion.data,
aes(x=inferred_meanDensity, y=baselinedOD, colour=Detector)) +
geom_point() +
scale_y_log10() +
scale_x_log10() +
labs(title = "data used for CFU to OD conversion")
# and again in terms of incubation temperature
ggplot(conversion.data,
aes(x=inferred_meanDensity, y=baselinedOD, colour=Incubation_temp)) +
geom_point() +
scale_y_log10() +
scale_x_log10() +
labs(title = "data used for CFU to OD conversion")
# a linear model with density as the predictor
CFU_to_OD.lm <- with(conversion.data,
lm(log10_OD ~ log10_density))
summary(CFU_to_OD.lm)##
## Call:
## lm(formula = log10_OD ~ log10_density)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.48332 -0.26051 -0.03599 0.25091 0.43137
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.28368 0.14465 -43.44 <2e-16 ***
## log10_density 0.89547 0.02383 37.58 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2701 on 336 degrees of freedom
## Multiple R-squared: 0.8078, Adjusted R-squared: 0.8073
## F-statistic: 1413 on 1 and 336 DF, p-value: < 2.2e-16# a function to convert from CFU to OD
convert_CFU_to_OD <- function(CFU) {
10^(coef(CFU_to_OD.lm)[1] + coef(CFU_to_OD.lm)[2]*log10(CFU))
}
# check what these predictions look like
some_CFU_vals <- seq(from=3.000e+00, to=4.362e+10, length.out=10)
pred_OD <- sapply(some_CFU_vals,
convert_CFU_to_OD)
ggplot(conversion.data,
aes(x=inferred_meanDensity, y=baselinedOD)) +
geom_point() +
scale_y_log10() +
scale_x_log10() +
geom_line(data=data.frame(some_CFU_vals, pred_OD), aes(x=some_CFU_vals, y=pred_OD)) +
labs(title = "checking CFU to OD conversion (line) against data (points)")
# clean up variables that will not be used in further analyses
rm(all_calib.data, conversion.data, pred_OD, some_CFU_vals)Analysis of CFU data
Sections I-II from protocols
High-throughput growth curves (exponential phase).docx and
High-throughput growth curves (stationary phase).docx where
used to create the single-use inoculum ‘shots’. Then, (within 2 months)
these shots were plated for CFU estimation by following the protocol
colony forming units (CFU) 6x6 drop plating.docx.
The CFU data is used to estimate the growth probability at different elevated temperatures (35C and 40C) as compared to normal lab culture conditions for soil bacteria (28C). It’s also used in conjunction with the growth curve data to inform us about the starting concentration of cells in each replicate for the time-to-detection method below, as indicated in Mytilinaios et al. 2015.
Load & parse the data
# load the raw data:
CFUraw_survival <- read_excel("./data/CFU_data_2023-12-13.xlsx", sheet="Temp_survival")
# ^^ this is the CFU's at different heat temperatures (35*C & 40*C) as compared to control (28*C). It will be analyzed in this section now.
CFUraw_inocula <- read_excel("./data/CFU_data_2023-12-13.xlsx", sheet="inocula")
# ^^ this is for the estimated inocula of the growth curves. It will be used in the analysis of growth curves below.
# annotate the Samples with their species and strain information
annot.df <- data.frame(Sample = c("BSC001", "BSC002", "BSC003", "BSC004", "BSC005", "BSC006", "BSC007", "BSC008", "BSC009", "BSC010", "BSC015", "BSC019", "CK101", "CK102", "CK103","CK104"),
Species = c(rep("P. putida", 3), rep("P. veronii", 3), rep("P. knackmussii", 2), "P. plecoglossicida", rep("P. putida", 2), "P. grimontii", rep("P. protegens", 4)),
Strain = c("F1", "KT2440", "uwc 2", rep(NA, 3), "B13", "B14", NA, "mt-2 KT2440", "KT2440", NA, rep("Pf5", 2), rep("CHAO", 2)),
Species_in_Expt2 = c(rep(TRUE, 3), rep(TRUE, 3), rep(FALSE, 2), FALSE, rep(TRUE, 2), TRUE, rep(TRUE, 4)))
CFUraw_survival <- inner_join(CFUraw_survival, annot.df, by=c("Sample"))
CFUraw_inocula <- inner_join(CFUraw_inocula, annot.df, by=c("Sample"))
# remove data points where the plate was noted to have dried out
CFUraw_inocula <- CFUraw_inocula[!with(CFUraw_inocula, grepl("Plate dried out after", Notes)),]
# remove data points where fewer than 15 cells were counted FOR THE 28*C CONTROL condition as these are not trustworthy estimates
CFUraw_survival <- CFUraw_survival %>% filter(CFU_at_28 > 14)
# calculate the survival under the extreme temperatures as a fraction of the density at 28
CFUraw_survival <- CFUraw_survival %>% mutate(survival = CFU_at_stress/CFU_at_28) %>%
rename(Temp = Stress_Temp) # rename this column to avoid confusion below
# get mean and bootstrapped confidence intervals
CFU_survival.df <- CFUraw_survival %>% group_by(Sample, Species, Strain, Temp, Inoculum, Species_in_Expt2) %>%
summarise(ci = list(mean_cl_boot(survival) %>%
rename(mean_surviv=y, lwr_ci=ymin, upr_ci=ymax))) %>% unnest(cols = c(ci))## `summarise()` has grouped output by 'Sample', 'Species', 'Strain', 'Temp',
## 'Inoculum'. You can override using the `.groups` argument.# clean up
rm(annot.df)Does stationary phase protect cells from heat?
Based on the classical microbiology literature, I would expect that stationary phase cultures would be more resistant to high temperatures as compared to exponential phase cultures. Let’s check if there’s a significant effect of the inoculum phase (i.e., exponential or stationary) on the CFU. To do this we will use a two-way ANOVA.
# convert the columns of the independent variables to factors
CFUraw_survival$Temp <- as.factor(CFUraw_survival$Temp)
CFUraw_survival$Species <- as.factor(CFUraw_survival$Species)
CFUraw_survival$Inoculum <- as.factor(CFUraw_survival$Inoculum)
# do the anova
# no interaction
survival.anova <- aov(survival ~ Inoculum + Species + Temp, data = CFUraw_survival)
# interaction
survival.anova.ALLinteractions <- aov(survival ~ Inoculum * Species * Temp, data = CFUraw_survival)
# use AIC to select the preferred model
model.set <- list(survival.anova, survival.anova.ALLinteractions)
model.names <- c("two.way", "interaction")
aictab(model.set, modnames = model.names)# display the results
summary(survival.anova.ALLinteractions)## Df Sum Sq Mean Sq F value Pr(>F)
## Inoculum 1 0.291 0.291 6.420 0.01222 *
## Species 5 18.113 3.623 79.869 < 2e-16 ***
## Temp 1 9.650 9.650 212.749 < 2e-16 ***
## Inoculum:Species 5 3.187 0.637 14.051 1.96e-11 ***
## Inoculum:Temp 1 0.015 0.015 0.329 0.56678
## Species:Temp 5 5.805 1.161 25.595 < 2e-16 ***
## Inoculum:Species:Temp 5 0.775 0.155 3.419 0.00576 **
## Residuals 164 7.439 0.045
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#TukeyHSD(survival.anova.ALLinteractions) # don't print this because it's very long and not so meaningful
# clean up
rm(survival.anova, survival.anova.ALLinteractions, model.set, model.names)So there’s a very slight effect of inoculum phase, but this is definitely less important than the effect of species. In fact, it seems that the inoculum effect actually interacts with the species. Neat!
Plot CFU data
# re-order the Species for appropriate colour plotting
CFU_survival.df$Species <- factor(CFU_survival.df$Species,
levels=c("P. veronii", "P. grimontii", "P. plecoglossicida", "P. protegens", "P. knackmussii", "P. putida"))
# for appropriate labelling of the facet with degree symbol, create a factorial version of Temp
CFU_survival.df$temp_as_factor <- factor(CFU_survival.df$Temp)
levels(CFU_survival.df$temp_as_factor) <- c("thirtyfive", "forty")
# for plotting little skulls to emphasize lack of growth
CFU_survival.df$Shape <- ifelse(CFU_survival.df$mean_surviv==0, "\u2620", "\u25cf")
plot.CFU_survivalEXP <- CFU_survival.df %>% filter(Inoculum == "Exponential") %>%
arrange(Temp, Species)
plot.CFU_survivalSTA <- CFU_survival.df %>% filter(Inoculum == "Stationary") %>%
arrange(Temp, Species)
print("exponential phase plot below:")## [1] "exponential phase plot below:"# plot the exponential phase data with Temp as rows
#ggplot(CFU_survival.df %>% filter(Inoculum == "Exponential"),
# aes(y=mean_surviv, x=Species, colour=Species)) +
# facet_grid(temp_as_factor ~ .,
# labeller = as_labeller(c(thirtyfive='35\u00b0C',
# forty='40\u00b0C'))) +
# geom_hline(yintercept = 1, colour="grey") +
# geom_point(position = position_dodge2(width = 0.8, preserve = "single"), alpha=0.5, shape=plot.CFU_survivalEXP$Shape, size=3) +
# geom_linerange(aes(ymin=lwr_ci, ymax=upr_ci), position = position_dodge2(width = 0.8, preserve = "single")) +
# scale_colour_manual(values=ALLsp_palette) +
# theme(axis.text.x=element_text(angle = 60, hjust = 0.95)) +
# labs(y="Effect Size (CFU stress/CFU 28*C)") +
# theme(strip.background = element_rect(fill = "white", colour="black"),
# strip.text = element_text(color = "black"))
# plot the exponential phase data again but with Temp as columns
ggplot(CFU_survival.df %>% filter(Inoculum == "Exponential"),
aes(y=mean_surviv, x=Species, colour=Species)) +
facet_wrap(. ~ temp_as_factor,
labeller = as_labeller(c(thirtyfive='35\u00b0C',
forty='40\u00b0C'))) +
geom_hline(yintercept = 1, colour="grey") +
geom_point(position = position_dodge2(width = 1), alpha=0.5, shape=plot.CFU_survivalEXP$Shape, size=3) +
geom_linerange(aes(ymin=lwr_ci, ymax=upr_ci), position = position_dodge2(width = 1)) +
scale_colour_manual(values=ALLsp_palette) +
theme(axis.text.x=element_text(angle = 60, hjust = 0.95)) +
labs(y="Heat Growth ± 95% CI (heat CFU / 28\u00b0C CFU)") +
theme(strip.background = element_rect(fill = "white", colour="black"),
strip.text = element_text(color = "black"),
#axis.title.y = element_text(size = 12),
axis.text.x=element_blank())## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_segment()`).
print("stationary phase plot below:")## [1] "stationary phase plot below:"# plot the stationary phase data with Temp as columns
ggplot(CFU_survival.df %>% filter(Inoculum == "Stationary"),
aes(y=mean_surviv, x=Species, colour=Species)) +
facet_wrap(. ~ temp_as_factor,
labeller = as_labeller(c(thirtyfive='35\u00b0C',
forty='40\u00b0C'))) +
geom_hline(yintercept = 1, colour="grey") +
geom_point(position = position_dodge2(width = 1), alpha=0.5, shape=plot.CFU_survivalSTA$Shape, size=3) +
geom_linerange(aes(ymin=lwr_ci, ymax=upr_ci), position = position_dodge2(width = 1)) +
scale_colour_manual(values=ALLsp_palette) +
theme(axis.text.x=element_text(angle = 60, hjust = 0.95)) +
labs(y="Heat Growth ± 95% CI (heat CFU / 28\u00b0C CFU)") +
theme(strip.background = element_rect(fill = "white", colour="black"),
strip.text = element_text(color = "black"),
#axis.title.y = element_text(size = 12),
axis.text.x=element_blank())
# clean up
rm(plot.CFU_survivalEXP, plot.CFU_survivalSTA, CFUraw_survival)Analysis of growth curves
Analysis of growth curve data that was generated by following the
protocols
High-throughput growth curves (exponential phase).docx and
High-throughput growth curves (stationary phase).docx.
By growing bacteria from a serial dilution, it is possible to estimate their intrinsic growth rate (\(\mu\)) by measuring the inoculum size and the time it takes for each dilution to cross some threshold (“time-to-detection”, abbreviated TTD), then rearranging the simple exponential growth equation \(N_t = N_0e^{\mu t}\). It’s also possible to see if there’s a lag time and to test if other models are a better fit to the data. See Mytilinaios et al. 2015. Remember that TTD is inversely related to growth rate: a fast growing strain will have a shorter TTD and a slow growing strain will have a longer TTD.
The calibration curve, CFUs, and heat growth data above are incorporated into the analysis.
Load & parse data
Load the data using the function get_ODdata_block
defined in “H1_extract_data_fun_v3-1.R”. Note that pairs of replicate
plates that were run on the same day (i.e., “blocks”) can be loaded
together in the same function call with the date as the block ID. The OD
data for each well is loaded. Also, the temperature measured at each
time-point by the instrument is given as the mean over the entire growth
curve rounded to the nearest integer.
Exclusion of contaminated wells: make sure that you exclude any contaminated wells from the data BEFORE exporting to text file. You will have to do this by noting which wells have unexpected fluorescence emission via fluorescent scanning of the plate with the BioTek H1 Synergy multimodal plate reader at the end of the growth curve.
#, eval=FALSE}
# this function runs with some warnings. Hermina looked over the resultant data.frame and it seems fine. Will fix this if they find some time...
# load the data
# list.files(pattern = 'C\\.txt$')
# Exponential Phase Inoculum
EX1raw_data.df <- get_ODdata_block(blockID = "2023-10-30",
filename_1 = "./data/2023-10-30_12strains_30C.txt",
filename_2 = "./data/2023-10-30_12strains_40C.txt")
EX2raw_data.df <- get_ODdata_block(blockID = "2023-11-01",
filename_1 = "./data/2023-11-01_12strains_25C.txt",
filename_2 = "./data/2023-11-01_12strains_35C.txt")
EX3raw_data.df <- get_ODdata_block(blockID = "2023-11-03",
filename_1 = "./data/2023-11-03_12strains_30C.txt",
filename_2 = "./data/2023-11-03_12strains_40C.txt")
EX4raw_data.df <- get_ODdata_block(blockID = "2023-11-06",
filename_1 = "./data/2023-11-06_12strains_25C.txt",
filename_2 = "./data/2023-11-06_12strains_35C.txt")
EX5raw_data.df <- get_ODdata_block(blockID = "2023-11-08",
filename_1 = "./data/2023-11-08_12strains_30C.txt")## [1] "Only loading 1 file for block 2023-11-08"# error with 2023-11-08_12strains_40C.txt
# when extract_OD_df splits the lines strings on tab and then tries to coerce them into a matrix, R complains that "data length [27354] is not a sub-multiple or multiple of the number of columns [98]"
# I can't find the issue and Excel can't either
# so just load this one file from csv
broken_txt <- read.csv("./data/2023-11-08_12strains_40C.csv")
# convert the wide data into long tidy data as usually done by get_OD_data
broken_txt$Time <- sapply(broken_txt$Time,
function(x) as.numeric(hms(x))/(60*60))
# now we can easily convert the whole data.frame to numeric
broken_txt[] <- lapply(broken_txt, as.numeric)
# convert Temp into the mean temperature across the run
broken_txt$Temp <- round(mean(broken_txt$Temp), digits=1)
# reshape the data to long format
broken_txt <- pivot_longer(broken_txt, cols=!c("Time", "Temp"), names_to="Well", values_to="OD")
# add block and plate columns
broken_txt$block <- "2023-11-08"
broken_txt$plate <- 2
# combine it with data from the rest of its corresponding block
# first need to change Temp factor back to numeric
EX5raw_data.df$Temp <- as.numeric(as.character(EX5raw_data.df$Temp))
EX5raw_data.df <- rbind(EX5raw_data.df, broken_txt)
# then change it back to a factor to prevent downstream issues
EX5raw_data.df$Temp <- as.factor(EX5raw_data.df$Temp)
rm(broken_txt)
EX6raw_data.df <- get_ODdata_block(blockID = "2023-11-10",
filename_1 = "./data/2023-11-10_12strains_25C.txt")## [1] "Only loading 1 file for block 2023-11-10" # there's no 2023-11-10_12strains_35C.txt because this file got corrupted
# Stationary Phase Inoculum
ST1raw_data.df <- get_ODdata_block(blockID = "2023-11-13",
filename_1 = "./data/2023-11-13_12strains_25C.txt",
filename_2 = "./data/2023-11-13_12strains_35C.txt")
ST2raw_data.df <- get_ODdata_block(blockID = "2023-11-20",
filename_1 = "./data/2023-11-20_12strains_30C.txt",
filename_2 = "./data/2023-11-20_12strains_40C.txt")
ST3raw_data.df <- get_ODdata_block(blockID = "2023-11-15",
filename_1 = "./data/2023-11-15_12strains_25C.txt",
filename_2 = "./data/2023-11-15_12strains_35C.txt")
ST4raw_data.df <- get_ODdata_block(blockID = "2023-11-22",
filename_1 = "./data/2023-11-22_12strains_30C.txt",
filename_2 = "./data/2023-11-22_12strains_40C.txt")
ST5raw_data.df <- get_ODdata_block(blockID = "2023-11-17",
filename_1 = "./data/2023-11-17_12strains_25C.txt",
filename_2 = "./data/2023-11-17_12strains_35C.txt")
ST6raw_data.df <- get_ODdata_block(blockID = "2023-11-24",
filename_1 = "./data/2023-11-24_12strains_30C.txt",
filename_2 = "./data/2023-11-24_12strains_40C.txt")Annotate the data using the function annotate_data
defined in “H1_extract_data_fun_v3-1.R”. The option
layout = "12columns" assumes the data is composed of 12
samples laid out over columns 1-12 and rows B-G are assumed to be a
dilution series from \(10^{-1}\) to
\(10^{-6}\) (rows A & H are assumed
to be blank).
#, eval=FALSE}
# define vectors for each of the 3 layouts
layout1 = c("Pseudomonas putida KT2440 - mScarlet","Pseudomonas putida uwc 2 - BFP","Pseudomonas putida F1 - sYFP","Pseudomonas veronii - mScarlet","Pseudomonas knackmussii B14 - GFP","Pseudomonas plecoglossicida - NA","Pseudomonas veronii - mCherry","Pseudomonas putida mt-2 KT2440 - NA","Pseudomonas veronii - BFP","Pseudomonas knackmussii B13 - mCherry","Pseudomonas putida KT2440 - GFP","Pseudomonas grimontii - NA")
layout2 = c("Pseudomonas veronii - mCherry","Pseudomonas veronii - BFP","Pseudomonas putida mt-2 KT2440 - NA","Pseudomonas knackmussii B13 - mCherry","Pseudomonas putida KT2440 - GFP","Pseudomonas grimontii - NA","Pseudomonas putida KT2440 - mScarlet","Pseudomonas putida F1 - sYFP","Pseudomonas putida uwc 2 - BFP","Pseudomonas veronii - mScarlet","Pseudomonas knackmussii B14 - GFP","Pseudomonas plecoglossicida - NA")
layout3 = c("Pseudomonas veronii - mScarlet","Pseudomonas knackmussii B14 - GFP","Pseudomonas plecoglossicida - NA","Pseudomonas putida KT2440 - mScarlet","Pseudomonas putida uwc 2 - BFP","Pseudomonas putida F1 - sYFP","Pseudomonas knackmussii B13 - mCherry","Pseudomonas grimontii - NA","Pseudomonas putida KT2440 - GFP","Pseudomonas veronii - mCherry","Pseudomonas veronii - BFP","Pseudomonas putida mt-2 KT2440 - NA")
# annotate the data for each raw data.frame
# Exponential Phase Inoculum
EX1raw_data.df <- annotate_data(EX1raw_data.df, layout = "12columns", samples=layout1)
EX2raw_data.df <- annotate_data(EX2raw_data.df, layout = "12columns", samples=layout1)
EX3raw_data.df <- annotate_data(EX3raw_data.df, layout = "12columns", samples=layout2)
EX4raw_data.df <- annotate_data(EX4raw_data.df, layout = "12columns", samples=layout2)
EX5raw_data.df <- annotate_data(EX5raw_data.df, layout = "12columns", samples=layout3)
EX6raw_data.df <- annotate_data(EX6raw_data.df, layout = "12columns", samples=layout3)
# Stationary Phase Inoculum
ST1raw_data.df <- annotate_data(ST1raw_data.df, layout = "12columns", samples=layout1)
ST2raw_data.df <- annotate_data(ST2raw_data.df, layout = "12columns", samples=layout1)
ST3raw_data.df <- annotate_data(ST3raw_data.df, layout = "12columns", samples=layout2)
ST4raw_data.df <- annotate_data(ST4raw_data.df, layout = "12columns", samples=layout2)
ST5raw_data.df <- annotate_data(ST5raw_data.df, layout = "12columns", samples=layout3)
ST6raw_data.df <- annotate_data(ST6raw_data.df, layout = "12columns", samples=layout3)
# combine exponential replicates into one data.frame
exponential.df <- rbind(EX1raw_data.df, EX2raw_data.df, EX3raw_data.df, EX4raw_data.df, EX5raw_data.df, EX6raw_data.df)
exponential.df$Inoculum <- "Exponential"
# combine stationary replicates into one data.frame
stationary.df <- rbind(ST1raw_data.df, ST2raw_data.df, ST3raw_data.df, ST4raw_data.df, ST5raw_data.df, ST6raw_data.df)
stationary.df$Inoculum <- "Stationary"
# combine all replicates into one giant data.frame
ALLraw_data.df <- rbind(exponential.df, stationary.df)
# remove the column titled Replicate as it is meaningless
ALLraw_data.df <- ALLraw_data.df %>% select(-Replicate)
# create a unique ID for each growth curve by pasting block, plate, and well
ALLraw_data.df$uniqID <- with(ALLraw_data.df, paste(Well, plate, block, Inoculum))
# the time column is giving very stupid and meaningless values because of division by 3
# round to nearest second
ALLraw_data.df$Time <- round(ALLraw_data.df$Time, digits = 4)
# the Gen5 program inserted ** around the values that Anjaney excluded as outliers.
# this leads to NA values that lead to downstream problems
# loop through all uniqID's to identify and remove replicates that are all NAs
for(curve in unique(ALLraw_data.df$uniqID)){
df <- ALLraw_data.df %>% filter(uniqID == curve)
if(all(is.na(df$OD))){ # identify replicates where all OD values are NA
ALLraw_data.df <-ALLraw_data.df %>% filter(uniqID != curve)
}
}
# cleanup
rm(df, curve)
# change all Temp==26 values into 25 because that's our treatment
ALLraw_data.df$Temp <- as.numeric(as.character(ALLraw_data.df$Temp))
ALLraw_data.df$Temp[which(ALLraw_data.df$Temp==26)] <- 25
# to keep your work space clean, you can now delete the many data.frames containing the individual blocks
rm(EX1raw_data.df, EX2raw_data.df, EX3raw_data.df, EX4raw_data.df, EX5raw_data.df, EX6raw_data.df, ST1raw_data.df, ST2raw_data.df, ST3raw_data.df, ST4raw_data.df, ST5raw_data.df, ST6raw_data.df, layout1, layout2, layout3, exponential.df, stationary.df)(Note that here plate 1 corresponds with H1 Synergy and plate 2 corresponds with Epoch.)
The P. protegens data was acquired by Hermina in Spring 2024, after Anjaney left. Load that data now:
#, eval=FALSE}
# load the data as separate data.frames
PpEXP1raw_data.df <- get_ODdata_block(blockID = "2024-03-13",
filename_1 = "./data/Pprotegens_exponential_40C--13March24.txt",
filename_2 = "./data/Pprotegens_exponential_35C--13March24.txt")
PpEXP2raw_data.df <- get_ODdata_block(blockID = "2024-03-11",
filename_1 = "./data/Pprotegens_exponential_30C--11March24.txt",
filename_2 = "./data/Pprotegens_exponential_25C--11March24.txt")
PpSTA1raw_data.df <- get_ODdata_block(blockID = "2024-03-15",
filename_1 = "./data/Pprotegens_stationary_40C--15March24.txt",
filename_2 = "./data/Pprotegens_stationary_30C--15March24.txt")
PpSTA2raw_data.df <- get_ODdata_block(blockID = "2024-03-18",
filename_1 = "./data/Pprotegens_stationary_35C--18March24.txt",
filename_2 = "./data/Pprotegens_stationary_25C--18March24.txt")
# annotate the data
layout1 <- c(rep("Pseudomonas protegens Pf5 - mTourquoise", 3), rep("Pseudomonas protegens Pf5 - mCherry", 3), rep("Pseudomonas protegens CHAO - GFP", 3), rep("Pseudomonas protegens CHAO - mCherry", 3))
layout2 <- c(rep(c("Pseudomonas protegens Pf5 - mCherry", "Pseudomonas protegens CHAO - GFP", "Pseudomonas protegens CHAO - mCherry", "Pseudomonas protegens Pf5 - mTourquoise"), times=2), "Pseudomonas protegens Pf5 - mCherry", "Pseudomonas protegens Pf5 - mTourquoise", "Pseudomonas protegens CHAO - mCherry", "Pseudomonas protegens CHAO - GFP")
layout3 <- rep(c("Pseudomonas protegens Pf5 - mTourquoise","Pseudomonas protegens Pf5 - mCherry","Pseudomonas protegens CHAO - GFP","Pseudomonas protegens CHAO - mCherry"), times=3)
# annotate the data for each raw data.frame
# Exponential Phase Inoculum
PpEXP1raw_data.df <- annotate_data(PpEXP1raw_data.df, layout = "12columns", samples=layout1)
PpEXP2raw_data.df <- annotate_data(PpEXP2raw_data.df, layout = "12columns", samples=layout2)
# Stationary Phase Inoculum
PpSTA1raw_data.df <- annotate_data(PpSTA1raw_data.df, layout = "12columns", samples=layout3)
PpSTA2raw_data.df <- annotate_data(PpSTA2raw_data.df, layout = "12columns", samples=layout2)
# combine exponential replicates into one data.frame
exponential.df <- rbind(PpEXP1raw_data.df, PpEXP2raw_data.df)
exponential.df$Inoculum <- "Exponential"
# combine stationary replicates into one data.frame
stationary.df <- rbind(PpSTA1raw_data.df, PpSTA2raw_data.df)
stationary.df$Inoculum <- "Stationary"
# combine all replicates into one giant data.frame
Ppraw_data.df <- rbind(exponential.df, stationary.df)
# remove the column titled Replicate as it is meaningless
Ppraw_data.df <- Ppraw_data.df %>% select(-Replicate)
# create a unique ID for each growth curve by pasting block, plate, and well
Ppraw_data.df$uniqID <- with(Ppraw_data.df, paste(Well, plate, block, Inoculum))
# the time column is giving very stupid and meaningless values because of division by 3
# round to nearest second
Ppraw_data.df$Time <- round(Ppraw_data.df$Time, digits = 4)
# the Gen5 program inserted ** around the values that Anjaney excluded as outliers.
# this leads to NA values that lead to downstream problems
# loop through all uniqID's to identify and remove replicates that are all NAs
for(curve in unique(Ppraw_data.df$uniqID)){
df <- Ppraw_data.df %>% filter(uniqID == curve)
if(all(is.na(df$OD))){ # identify replicates where all OD values are NA
Ppraw_data.df <- Ppraw_data.df %>% filter(uniqID != curve)
}
}
# cleanup
rm(df, curve)
# change all Temp==26 values into 25 because that's our treatment
Ppraw_data.df$Temp <- as.numeric(as.character(Ppraw_data.df$Temp))
Ppraw_data.df$Temp[which(Ppraw_data.df$Temp==26)] <- 25
# combine the P. protegens data with the rest of the data...
ALLraw_data.df <- rbind(ALLraw_data.df, Ppraw_data.df)
# clean up
rm(PpEXP1raw_data.df, PpEXP2raw_data.df, PpSTA1raw_data.df, PpSTA2raw_data.df, annotate_data, extract_OD_df, get_OD_data, get_ODdata_block, layout1, layout2, layout3, exponential.df, stationary.df, Ppraw_data.df)Do the baseline subtraction to get the OD values: baseline subtract using the median OD value of the blanks for each time point, for each plate.
Then let’s plot the baseline subtracted data OD for all culture phases, samples, and temperatures. The colours show growth curves with different dilutions, from \(10^{-1}\) to \(10^{-6}\) (note that these come from different batches, i.e., that were grown on different days). The black horizontal line shows the detection threshold (baselined OD = 0.05) that is used for the TTD analysis below.
#, eval=FALSE}
# baseline subtract using the median OD value of the blanks for each time point, for each plate
# create a data.frame containing only the blank samples for all plates
data.blanks <- ALLraw_data.df %>% filter(Sample=="BLK")
# calculate the median OD value for the blanks at each time point for each plate.
medianOD <- data.blanks %>% group_by(Time, block, plate) %>% # group the blanks together for each time point and plate
summarise(medianBlankOD = median(OD, na.rm=TRUE)) # get median OD value## `summarise()` has grouped output by 'Time', 'block'. You can override using the
## `.groups` argument. # put the data.blanks information back into the full dataset
ALL_data.df <- inner_join(ALLraw_data.df %>% filter(Sample!="BLK"), # remove blank samples
medianOD)## Joining with `by = join_by(Time, block, plate)` # Finally, do the baseline subtraction by subtracting the medianBlankOD from the OD values
ALL_data.df$baselinedOD <- ALL_data.df$OD - ALL_data.df$medianBlankOD
# plot the data
ggplot(filter(ALL_data.df, Sample == "Pseudomonas protegens Pf5 - mTourquoise", Inoculum == "Stationary", Temp == "25"),
aes(y= baselinedOD, x=Time, colour=Dilution, group=uniqID )) +
scale_y_log10() +
geom_line() +
geom_hline(yintercept = 0.05) +
theme(legend.position = "none") +
labs(title="P. protegens, Stationary, 25C, TTDCutoff = 0.05", x="Time (hrs)", y="A600 (corrected O.D.)") ## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_line()`).
for(phase in unique(ALL_data.df$Inoculum)){
for(smpls in levels(ALL_data.df$Sample)[levels(ALL_data.df$Sample) != "BLK"]){
for(temps in unique(ALL_data.df$Temp)) {
temp_data_for_plotting <- ALL_data.df %>%
filter(Inoculum == phase,
Sample == smpls,
Temp == temps)
plot <- ggplot(temp_data_for_plotting,
aes(y = baselinedOD, x = Time, colour = Dilution, group = uniqID)) +
scale_y_log10() +
geom_line() +
geom_hline(yintercept = 0.05) +
theme(legend.position = "none") +
labs(title = paste0(phase, " at ", temps, "*C, ", smpls),
x="Time (hrs)", y="A600 (corrected O.D.)")
print(plot)
}
}
}## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 382 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 2766 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 303 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 157 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1208 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1275 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 845 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 230 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 211 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1071 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 339 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 453 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 497 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 2598 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 194 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 117 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 672 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1208 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 621 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 285 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 592 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 662 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 141 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 113 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 154 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 817 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 213 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 119 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 214 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 704 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 84 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 124 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 256 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 565 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 374 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 178 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 238 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 3194 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 154 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 513 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 200 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 2490 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 557 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 862 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 160 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1812 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 637 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1782 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 306 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 872 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 571 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 759 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 297 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 2602 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 918 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 560 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 356 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1867 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 545 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 582 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 519 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 463 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 456 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 339 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 269 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1794 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 327 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 473 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 951 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1198 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 783 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1278 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 238 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 174 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 105 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 101 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1568 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 163 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 301 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 319 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 863 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 446 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 89 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 339 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 587 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 178 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 352 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 158 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 654 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 401 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 162 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 303 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 198 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 286 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 263 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 118 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 845 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 216 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 146 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 678 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 2597 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 606 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1869 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 494 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 2179 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 485 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 2556 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 357 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1836 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 446 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1430 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 368 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 748 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 381 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 233 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 146 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1322 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 195 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 126 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 418 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 1035 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 153 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 197 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 325 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 353 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 37 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 208 rows containing missing values or values outside the scale range
## (`geom_line()`).
# clean-up
rm(data.blanks, medianOD, plot, temp_data_for_plotting, phase, smpls, temps)Improve on the annotation of the Samples to create dedicated columns for species, strain, and fluorescence and also to add the sample ID in the TERR collection.
Then, let’s see how many data points we have for all samples and treatments…
Note that due to some error with the plate reader, the 3rd replicate of 35*C (in block 2023-11-10) was corrupted so we’re missing those data points :(
#, eval=FALSE}
# parse the details stored in the Sample column into dedicated columns...
ALL_data.df$temp_Sample <- as.character(ALL_data.df$Sample)
# separate the species, strain, and fluorophore information from the unique sample names
ALL_data.df <- ALL_data.df %>% separate(col=temp_Sample, sep=" - ", into=c("Genotype", "Fluor"))
# parse the Genotype column into a species column
ALL_data.df$Species <- sub("Pseudomonas (\\w+)", "\\1",
str_extract(ALL_data.df$Genotype, "Pseudomonas (\\w+)"))
ALL_data.df$Species <- paste("P.", ALL_data.df$Species)
# For some samples, we have additional strain information. Parse this from the Genotype column into a strain column
ALL_data.df$Strain <- NA
ALL_data.df$Strain[as.logical(str_count(ALL_data.df$Genotype, "\\w+\\s\\w+\\s.*"))] <- sub("(\\S+\\s+){2}(.*)", "\\2", ALL_data.df$Genotype[as.logical(str_count(ALL_data.df$Genotype, "\\w+\\s\\w+\\s.*"))])
# Replace the information in the Sample column with the unique sample ID from our strain collection
ALL_data.df$Sample <- as.character(ALL_data.df$Sample)
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas putida F1 - sYFP"] <- "BSC001"
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas putida KT2440 - mScarlet"] <- "BSC002"
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas putida uwc 2 - BFP"] <- "BSC003"
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas veronii - BFP"] <- "BSC004"
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas veronii - mScarlet"] <- "BSC005"
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas veronii - mCherry"] <- "BSC006"
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas knackmussii B13 - mCherry"] <- "BSC007"
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas knackmussii B14 - GFP"] <- "BSC008"
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas plecoglossicida - NA"] <- "BSC009"
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas putida mt-2 KT2440 - NA"] <- "BSC010"
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas putida KT2440 - GFP"] <- "BSC015"
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas grimontii - NA"] <- "BSC019"
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas protegens Pf5 - mTourquoise"] <- "CK101"
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas protegens Pf5 - mCherry"] <- "CK102"
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas protegens CHAO - GFP"] <- "CK103"
ALL_data.df$Sample[ALL_data.df$Sample == "Pseudomonas protegens CHAO - mCherry"] <- "CK104"
# remove the genotype column
ALL_data.df <- ALL_data.df %>% select(-Genotype)
# check that all the data we loaded is actually there
check_data <- ALL_data.df %>% filter(Time < 0.2) %>%
group_by(Sample, Temp, Dilution, Inoculum) %>%
summarise(Count=n()) %>%
arrange(Inoculum, Sample, Temp, desc(Dilution)) %>% as.data.frame()## `summarise()` has grouped output by 'Sample', 'Temp', 'Dilution'. You can
## override using the `.groups` argument.# let's take a peak at the df...
check_data# okay that table is waaaaay too long, let's see an overview of it:
summary(check_data)## Sample Temp Dilution Inoculum
## Length:768 Min. :25.00 1e-06:128 Length:768
## Class :character 1st Qu.:28.75 1e-05:128 Class :character
## Mode :character Median :32.50 1e-04:128 Mode :character
## Mean :32.50 0.001:128
## 3rd Qu.:36.25 0.01 :128
## Max. :40.00 0.1 :128
## Count
## Min. :2.000
## 1st Qu.:3.000
## Median :3.000
## Mean :2.891
## 3rd Qu.:3.000
## Max. :3.000# 35*C in block 2023-11-10 was corrupted so that's why we're missing those data points
# Display all data that have fewer than 3 replicates
check_data %>% filter(Count < 3)# cleanup
rm(check_data, ALLraw_data.df)Data processing: Approximate the TTD
See from the data plots above (and Anjaney’s presentation) that a baselined OD value of 0.05 is a reasonable cut-off for the time-to-detection (TTD) such that all curves are still in exponential growth at this value (i.e., growth is linear when plotted on a log scale). However, this low value doesn’t work properly for wells whose baseline from the start of the experiment just happened to be higher than 0.05. That’s why I’m also setting a high TTD value of 0.13 for comparison.
We will use kernel smoothing to get a more fine-grained approximation for the TTD and find the exact time when the threshold is passed (i.e., as compared to Anjaney’s presentation, where he got a coarse-grained estimate of the TTD by finding the first sampled time point when baselined OD > 0.05).
#, eval=FALSE}
# set the cut-off to 0.05 because this is the value that Anjaney found where all growth curves are still in the "log-linear" phase.
TTDcutoff.low <- 0.05
tempdf <- ALL_data.df %>% filter(uniqID == "B1 1 2023-10-30 Exponential")
# perform linear interpolation using OD as the independent variable and Time as the dependent
tempdf.Time <- approx(y=tempdf$Time, x=tempdf$baselinedOD,
xout=TTDcutoff.low, ties="mean")$y #when there are non-unique values for baselinedOD, it takes the mean across these (that's why I'm NOT using the log of baselinedOD because then ties would be resolved using the geometric mean and I don't want that)
# plot to check that this is working properly
ggplot(tempdf,
aes(y= baselinedOD, x=Time)) +
scale_y_log10() +
geom_line() +
xlim(0,5) +
geom_hline(yintercept = TTDcutoff.low, col="blue") +
geom_vline(xintercept = tempdf.Time, col="red") +
theme(legend.position = "none") +
labs(title="Sanity check of TTD loess (kernal-smoothing)", x="Time (hrs)", y="A600 (corrected O.D.)") ## Warning: Removed 254 rows containing missing values or values outside the scale range
## (`geom_line()`).
rm(tempdf, tempdf.Time)
# initialize an empty data.frame for storage
TTD.df <- ALL_data.df %>% select(-Time, -Well, -OD, -medianBlankOD, -baselinedOD) %>% distinct()
TTD.df$TTD <- NA
# loop through all uniqID's
for(curve in unique(ALL_data.df$uniqID)) {
# calculate the detection time using interpolation
detect_t <- with(ALL_data.df %>% filter(uniqID == curve),
approx(y=Time, x=baselinedOD, xout=TTDcutoff.low, ties="mean")$y)
# save the time to the storage data.frame
TTD.df$TTD[TTD.df$uniqID == curve] <- detect_t
}
rm(curve, detect_t)
# check that NA values for TTD correspond with curves that failed to grow
# keep OD values for the unique ID's that have TTD equal to NA
tempdf <- ALL_data.df %>% filter(uniqID %in% TTD.df$uniqID[is.na(TTD.df$TTD)])
ggplot(tempdf,
aes(y = baselinedOD, x = Time, colour = uniqID, group = uniqID)) +
scale_y_log10() +
geom_line() +
#xlim(0,5) +
geom_hline(yintercept = TTDcutoff.low, col="blue") +
theme(legend.position = "none") +
labs(title="Checking growth curves with TTD 0.05=NA", x="Time (hrs)", y="A600 (corrected O.D.)")## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 51457 rows containing missing values or values outside the scale range
## (`geom_line()`).
ggplot(tempdf[which(tempdf$baselinedOD > TTDcutoff.low),],
aes(y = baselinedOD, x = Time, colour = uniqID, group = uniqID)) +
scale_y_log10() +
geom_line() +
geom_hline(yintercept = TTDcutoff.low, col="blue") +
geom_hline(yintercept = 0.13, col="red") +
labs(title="Check NA growth curves with TTD threshold=0.05", x="Time (hrs)", y="A600 (corrected O.D.)")
# let's repeat the TTD value calculation for slightly higher TTD cutoff where all growth curves can be included
# rename the first TTD value calculation in the storage data.frame
TTD.df <- TTD.df %>% rename(TTD_0.05 = TTD)
# define the higher TTD cutoff value and initialize storage vector
TTDcutoff.hi <- 0.13
TTD.df$TTD_0.13 <- NA
# loop through all uniqID's
for(curve in unique(ALL_data.df$uniqID)) {
# calculate the detection time using interpolation
detect_t <- with(ALL_data.df %>% filter(uniqID == curve),
approx(y=Time, x=baselinedOD, xout=TTDcutoff.hi, ties="mean")$y)
# save the time to the storage data.frame
TTD.df$TTD_0.13[TTD.df$uniqID == curve] <- detect_t
}
# check again that the NA TTD values represent growth curves that failed to grow
tempdf <- ALL_data.df %>% filter(uniqID %in% TTD.df$uniqID[is.na(TTD.df$TTD_0.13)])
ggplot(tempdf,
aes(y = baselinedOD, x = Time, colour = uniqID, group = uniqID)) +
scale_y_log10() +
geom_line() +
#xlim(0,5) +
geom_hline(yintercept = TTDcutoff.hi, col="red") +
theme(legend.position = "none") +
labs(title="Check NA growth curves with TTD threshold=0.13", x="Time (hrs)", y="A600 (corrected O.D.)")## Warning in transformation$transform(x): NaNs produced## Warning in scale_y_log10(): log-10 transformation introduced infinite values.## Warning: Removed 52225 rows containing missing values or values outside the scale range
## (`geom_line()`).
# check for agreement in the growth rates estimated by the 2 thresholds
ggplot(TTD.df, aes(x=TTD_0.05, y=TTD_0.13)) +
geom_abline(slope=1, intercept=0, colour="grey") +
geom_point(alpha=0.1) +
geom_smooth(se=FALSE, method="lm") +
labs(x="TTD for threshold = 0.05", y="TTD threshold = 0.13",
title="grey line: y=x")## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 482 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 482 rows containing missing values or values outside the scale range
## (`geom_point()`).
# display the correlation:
summary(lm(TTD.df$TTD_0.13 ~ TTD.df$TTD_0.05))##
## Call:
## lm(formula = TTD.df$TTD_0.13 ~ TTD.df$TTD_0.05)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3278 -0.2861 -0.0842 0.1609 12.2079
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.602497 0.029144 20.67 <2e-16 ***
## TTD.df$TTD_0.05 1.038598 0.002043 508.30 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5878 on 1736 degrees of freedom
## (482 observations deleted due to missingness)
## Multiple R-squared: 0.9933, Adjusted R-squared: 0.9933
## F-statistic: 2.584e+05 on 1 and 1736 DF, p-value: < 2.2e-16# save this intermediate data to output for reproducibility & easy loading in downstream analyses
save(ALL_data.df, TTD.df,
file="all_growthcurve_data.RData")
#clean up
rm(curve, detect_t, tempdf, TTDcutoff.hi, TTDcutoff.low)It seems that the lower TTD value (0.05) is good for most of the growth curves but there are a couple of growth curves that start with too high OD, even after baselining. That’s why I tried setting a higher TTD value (0.13). The problem with this higher value is that we lose some of the growth curves that only cross this threshold later in growth. shrug This is not a perfect method… But hopefully we can make some reasonable interpretations.
Below we can see that although bigger values are systematically estimated for TTD when the cut-off is bigger, there’s actually very good correlation between the estimates. In general I will prefer the lower cut-off value (0.05).
Before we combine the TTD data with the inoculum sizes from the CFU data, let’s convince ourselves that these two types of data are generally in agreement with one another. Use the data from both the higher and the lower TTD cutoffs as a binary indicator of growth (i.e., if we failed to find you at both cutoffs, then we assume that you didn’t grow at all). AKA survival at different temperatures. Then we will compare the binary indicator of growth with whether growth was observed or not by CFU.
#load("all_growthcurve_data.RData")
binary_survival <- TTD.df %>% filter(Temp > 25)
# convert the data TTD data into binary survival for both thresholds
binary_survival$TTD_0.05 <- !is.na(binary_survival$TTD_0.05)
binary_survival$TTD_0.13 <- !is.na(binary_survival$TTD_0.13)
# let's condense the TTD binary data by assuming that if you returned TRUE in either of the thresholds, then you grew.
# this is more conservative than strictly requiring growth in both TTD's...
binary_survival$TTD <- binary_survival$TTD_0.05 | binary_survival$TTD_0.13
# summarize the binary growth data into fraction of all replicates when growth was observed
binary_survival <- binary_survival %>% select(Sample, Species, Temp, Inoculum, Dilution, TTD) %>% # keep just the meta-data and the TTD binary survival
group_by(Sample, Temp, Inoculum) %>%
summarise(fraction_batch=mean(TTD))## `summarise()` has grouped output by 'Sample', 'Temp'. You can override using
## the `.groups` argument.# combine the fraction of growing batch cultures with the CFU data
temp_CFUsurviv <- CFU_survival.df %>% ungroup() %>% select(Sample, Temp, Inoculum, mean_surviv)
temp_CFUsurviv <- rbind(temp_CFUsurviv,
temp_CFUsurviv[temp_CFUsurviv$Temp == 35, ] %>% mutate(Temp = 30, mean_surviv = 1))
binary_survival <- inner_join(binary_survival,
temp_CFUsurviv,
by = c("Temp", "Sample", "Inoculum"))
# CFU growth > 1 is not particularly meaningful for this analysis
binary_survival$mean_surviv[binary_survival$mean_surviv > 1] <- 1
ggplot(binary_survival,
aes(y=mean_surviv, x=fraction_batch)) +
geom_point(alpha=0.1) +
labs(y="Fraction of CFUs that grew",
x="Fraction of batch cultures that grew")
# display the correlation:
summary(lm(binary_survival$mean_surviv ~ binary_survival$fraction_batch))##
## Call:
## lm(formula = binary_survival$mean_surviv ~ binary_survival$fraction_batch)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.77894 -0.05907 0.02093 0.06036 0.54028
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.05963 0.03826 -1.559 0.122
## binary_survival$fraction_batch 1.03870 0.04609 22.535 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1639 on 94 degrees of freedom
## Multiple R-squared: 0.8438, Adjusted R-squared: 0.8421
## F-statistic: 507.8 on 1 and 94 DF, p-value: < 2.2e-16# clean up
rm(temp_CFUsurviv, binary_survival)Batch culture and CFU data are not exactly the same but they tend to be in good agreement (\(R^2 \approx 0.84\)).
Next we will combine the estimated inoculum sizes with the TTD data in order to infer the growth rate.
# calculate the concentration of the cryostocked inoculum in the absence of stress (i.e, 28*C)
inoc_data <- CFUraw_inocula %>% filter(Temp == 28) %>% select(-Temp) %>%
mutate(Concentration = CFU/Vol_Plated_mL * 1/Dilution_Counted)
# get mean and bootstrapped confidence intervals for the cryostocked inoculum density
inoc_data <- inoc_data %>% group_by(Sample, Inoculum, Species, Strain, Species_in_Expt2) %>%
summarise(ci = list(mean_cl_boot(Concentration) %>%
rename(mean_inoc_den=y, lwrCI_inoc_den=ymin, uprCI_inoc_den=ymax))) %>% unnest(cols = c(ci))## `summarise()` has grouped output by 'Sample', 'Inoculum', 'Species', 'Strain'.
## You can override using the `.groups` argument.# now use the mean survival rates from CFU_survival.df to estimate the inoculum density at different temperatures
temp_CFUsurviv <- CFU_survival.df %>% select(-temp_as_factor, -Shape)
# assume that 25*C and 30*C have the same cell counts as 28*C ... in other words, survival = 1
temp_CFUsurviv$Temp <- temp_CFUsurviv$Temp - 10
temp_CFUsurviv$mean_surviv <- temp_CFUsurviv$lwr_ci <- temp_CFUsurviv$upr_ci <- 1
temp_CFUsurviv <- rbind(CFU_survival.df %>% select(-temp_as_factor, -Shape),
temp_CFUsurviv)
# combine the CFU survival data with the inoculum density data
inoc_data <- inner_join(inoc_data, temp_CFUsurviv,
by=c("Sample", "Species", "Strain", "Inoculum", "Species_in_Expt2"))
# finally, multiply the inoculum density with the CFU survival
inoc_data <- inoc_data %>%
mutate(mean_inoc_den = mean_inoc_den * mean_surviv,
lwrCI_inoc_den = lwrCI_inoc_den * lwr_ci,
uprCI_inoc_den = uprCI_inoc_den * upr_ci) %>%
select(-mean_surviv, -lwr_ci, -upr_ci)
# infer the inoculum density estimates for the dilution series
inoc_data <- do.call("rbind", replicate(6, inoc_data, simplify = FALSE))
inoc_data$Dilution <- rep(10^(-1*(1:6)), each=nrow(inoc_data)/6)
# calculate the mean N0
inoc_data <- inoc_data %>%
mutate(mean_N0 = mean_inoc_den * Dilution,
lwrCI_N0 = lwrCI_inoc_den * Dilution,
uprCI_N0 = uprCI_inoc_den * Dilution) %>%
select(-mean_inoc_den, -lwrCI_inoc_den, -uprCI_inoc_den)
# combine the data sets
# convert columns to compatible classes between the 2 data frames
TTD.df$Dilution <- as.numeric(levels(TTD.df$Dilution))[TTD.df$Dilution]
TTD.df <- inner_join(TTD.df, inoc_data,
by=c("Temp", "Sample", "Inoculum", "Species", "Strain", "Dilution"))
# exclude growth curves that started with fewer than 20 cells (too much stochasticity)
TTD.df <- TTD.df %>% filter(mean_N0 > 20)
###############################
# convert TTD into growth rate estimates
###############################
# estimate ln(N_t / N_0) in units of OD
TTD.df$inferredOD_meanN0 <- sapply(TTD.df$mean_N0,
convert_CFU_to_OD)
TTD.df <- TTD.df %>% mutate(ln_Ntsmall_over_N0 = log(0.05 / inferredOD_meanN0),
ln_Ntbig_over_N0 = log(0.13 / inferredOD_meanN0))
ggplot(TTD.df %>% filter(Inoculum=="Exponential"),
aes(x = ln_Ntsmall_over_N0, y = TTD_0.05, group=Temp, colour=as.factor(Temp))) +
facet_wrap(vars(Sample)) +
geom_point(alpha=0.5) +
geom_smooth(se=FALSE, method="lm", size=0.2) +
labs(title="Exponential; TTD cut-off 0.05",
x="ln(0.05 / OD converted inoculum)",
y="Time to Detection (hrs)",
colour="*C")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 17 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 17 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggplot(TTD.df %>% filter(Inoculum=="Exponential"),
aes(x=ln_Ntbig_over_N0, y=TTD_0.13, group=Temp, colour=as.factor(Temp))) +
facet_wrap(vars(Sample)) +
geom_point(alpha=0.5) +
geom_smooth(se=FALSE, method="lm", size=0.2) +
labs(title="Exponential; TTD cut-off 0.13",
x="ln(0.13 / OD converted inoculum)",
y="Time to Detection (hrs)",
colour="*C")## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 19 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 19 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggplot(TTD.df %>% filter(Inoculum=="Stationary"),
aes(x=ln_Ntsmall_over_N0, y=TTD_0.05, group=Temp, colour=as.factor(Temp))) +
facet_wrap(vars(Sample)) +
geom_point(alpha=0.5) +
geom_smooth(se=FALSE, method="lm", size=0.2) +
labs(title="Stationary; TTD cut-off 0.05",
x="ln(0.05 / OD converted inoculum)",
y="Time to Detection (hrs)",
colour="*C")## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 29 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 29 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggplot(TTD.df %>% filter(Inoculum=="Stationary"),
aes(x=ln_Ntbig_over_N0, y=TTD_0.13, group=Temp, colour=as.factor(Temp))) +
facet_wrap(vars(Sample)) +
geom_point(alpha=0.5) +
geom_smooth(se=FALSE, method="lm", size=0.2) +
labs(title="Stationary; TTD cut-off 0.13",
x="ln(0.13 / OD converted inoculum)",
y="Time to Detection (hrs)",
colour="*C")## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 18 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 18 rows containing missing values or values outside the scale range
## (`geom_point()`).
## EVERYTHING LOOKS GOOD SO NOW WE CAN ACTUALLY FIT THE LM'S AND EXTRACT SUMMARY STATISTICS FROM THEM:
# initialize variable for storage
Dil_growthrates.df <- data.frame(matrix(nrow=0, ncol=11))
# calculate the slope, sd, and r-squared of the slope
for(sam in unique(TTD.df$Sample)){
for(t in unique(TTD.df$Temp)){
for(inoc in unique(TTD.df$Inoculum)){
temp.df <- TTD.df %>% filter(Sample==sam, Temp==t, Inoculum==inoc)
if(sum(temp.df$mean_N0, na.rm=TRUE)>0 # if CFU data indicates cells can survive
& # and
sum(!is.na(temp.df$TTD_0.05))>4){ # there are 5 or more non-NA TTD values
# fit the linear model on ln(N_t / N_0) for N0 converted from CFU to OD
tmpOD_lm <- lm(TTD_0.05 ~ ln_Ntsmall_over_N0, temp.df)
# bootstrap the data to get the 95% CI of the slope
tmpOD_boot <- bootSlope(TTD_0.05 ~ ln_Ntsmall_over_N0, data = temp.df,
B = 1000, plot.hist = FALSE)
# estimate biomass growth rate from WITHOUT taking into account inocula CFU
tmp_lm <- lm(TTD_0.05 ~ log10(Dilution), temp.df)
# bootstrap the data to get the 95% CI of the slope
tmp_boot <- bootSlope(TTD_0.05 ~ log10(Dilution), data = temp.df,
B = 1000, plot.hist = FALSE)
# store the values
Dil_growthrates.df <- rbind(Dil_growthrates.df,
c(sam, t, inoc,
# summary statistics for the lm with converted OD
unname(1/coef(tmpOD_lm)[2]), # mean growth rate
unname(1/(quantile(tmpOD_boot, probs = 0.975,
na.rm = TRUE))), # lower 95% CI of growth rate
unname(1/(quantile(tmpOD_boot, probs = 0.025,
na.rm = TRUE))), # upper 95% CI of growth rate
summary(tmpOD_lm)$r.squared, # multiple R-squared
# also get the summary statistics for the lm on raw CFU
unname(exp(coef(tmp_lm)[2])), # mean growth rate
unname(exp(quantile(tmp_boot, probs = 0.025,
na.rm = TRUE))), # lower 95% CI of growth rate
unname(exp(quantile(tmp_boot, probs = 0.975,
na.rm = TRUE))), # upper 95% CI of growth rate
summary(tmp_lm)$r.squared) # multiple R-squared
)
}
if(sum(temp.df$mean_N0, na.rm=TRUE)==0 # if CFU data indicates cells CANNOT survive
| # or
sum(!is.na(temp.df$TTD_0.05))<5){ # if there are less than 5 non-NA TTD values
# store 0 values
Dil_growthrates.df <- rbind(Dil_growthrates.df,
c(sam, t, inoc, rep(0, 8))
)
# initialize variables to avoid unnecessary warnings
tmp_lm <- tmp_boot <- tmpOD_lm <- tmpOD_boot <- NA
}
rm(tmp_lm, tmp_boot, temp.df, tmpOD_lm, tmpOD_boot)
}
}
}##
## 1 observation(s) removed due to missing values.
##
## 1 observation(s) removed due to missing values.
##
## 3 observation(s) removed due to missing values.## Warning in cor(newx, newy): the standard deviation is zero## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero##
## 3 observation(s) removed due to missing values.## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero##
## 4 observation(s) removed due to missing values.
##
## 4 observation(s) removed due to missing values.## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero##
## 2 observation(s) removed due to missing values.
##
## 2 observation(s) removed due to missing values.
##
## 3 observation(s) removed due to missing values.
##
## 3 observation(s) removed due to missing values.## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero##
## 3 observation(s) removed due to missing values.
##
## 3 observation(s) removed due to missing values.
##
## 4 observation(s) removed due to missing values.## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero##
## 4 observation(s) removed due to missing values.## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero##
## 3 observation(s) removed due to missing values.## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero##
## 3 observation(s) removed due to missing values.## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero
## Warning in cor(newx, newy): the standard deviation is zero##
## 3 observation(s) removed due to missing values.
##
## 3 observation(s) removed due to missing values.## Warning in cor(newx, newy): the standard deviation is zero##
## 1 observation(s) removed due to missing values.
##
## 1 observation(s) removed due to missing values.# update the column names
colnames(Dil_growthrates.df) <- c("Sample", "Temp", "Inoculum",
"mu", "mu_lwrCI", "mu_uprCI", "r_sq",
"biomass_mu", "biomass_mu_lwrCI", "biomass_mu_uprCI", "biomass_rsq")
# change Temp to a numeric from character
Dil_growthrates.df$Temp <- as.numeric(Dil_growthrates.df$Temp)
# change mu to a numeric from character
Dil_growthrates.df$mu <- as.numeric(Dil_growthrates.df$mu)
Dil_growthrates.df$biomass_mu <- as.numeric(Dil_growthrates.df$biomass_mu)
# change mu_lwrCI to a numeric from character
Dil_growthrates.df$mu_lwrCI <- as.numeric(Dil_growthrates.df$mu_lwrCI)
Dil_growthrates.df$biomass_mu_lwrCI <- as.numeric(Dil_growthrates.df$biomass_mu_lwrCI)
# change mu_uprCI to a numeric from character
Dil_growthrates.df$mu_uprCI <- as.numeric(Dil_growthrates.df$mu_uprCI)
Dil_growthrates.df$biomass_mu_uprCI <- as.numeric(Dil_growthrates.df$biomass_mu_uprCI)
# change r_sq to a numeric from character
Dil_growthrates.df$r_sq <- as.numeric(Dil_growthrates.df$r_sq)
Dil_growthrates.df$biomass_rsq <- as.numeric(Dil_growthrates.df$biomass_rsq)
# add species/strain annotation
Dil_growthrates.df <- left_join(Dil_growthrates.df,
TTD.df %>% select(Sample, Species, Strain, Fluor, Species_in_Expt2) %>% distinct())## Joining with `by = join_by(Sample)`# check the r-squared values
with(Dil_growthrates.df, hist(r_sq))
# batch cultures that failed to grow reliably (need to have n=5 at minimum, corresponding with growth in at least 28% or 42% of incoulated wells) have 0 values for slope and r-squared.
# plot just the linear regressions that had low r-squared (there are 9 of them in total)
temp.low_rsq <- Dil_growthrates.df %>% filter(r_sq < 0.9) %>% filter(r_sq > 0) %>% select(Sample, Temp, Inoculum)
temp.low_rsq <- inner_join(temp.low_rsq, TTD.df) %>%
unite("sam_temp_inoc", c(Sample, Temp, Inoculum), remove=FALSE)## Joining with `by = join_by(Sample, Temp, Inoculum)` # plot
ggplot(temp.low_rsq,
aes(x=ln_Ntsmall_over_N0, y=TTD_0.05)) +
facet_wrap(vars(sam_temp_inoc)) +
geom_point(alpha=0.5) +
geom_smooth(se=FALSE, method="lm", size=0.2) +
labs(title="lm: r-sq < 0.85",
x="ln(0.05 / OD converted inoculum)",
y="Time to Detection (hrs)")## `geom_smooth()` using formula = 'y ~ x'## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_smooth()`).## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
# unfortunately, the r-squared is not a robust statistic of which linear regressions look unacceptable...
# clean up
rm(temp_CFUsurviv, CFU_survival.df, CFUraw_inocula, inoc_data, inoc, sam, t, temp.low_rsq)I removed all growth curves whose inoculated population size (by CFU estimates) was less than 20 cells. This is because low inocula will lead to very stochastic time-to-detection values, and we aren’t interested in that.
To get the 95% CI on the slope, I bootstrapped the data. Therefore I’m not very worried if there’s a weird outlier at the far end of the linear regression. This will ofc impact the mean estimate for the slope, but should be very robust to the CI.
Plot & Analyze
# for plotting little skulls to emphasize lack of growth
Dil_growthrates.df$Shape <- ifelse(Dil_growthrates.df$mu==0, "\u2620", "\u25cf")
# re-organize the data.frame for correct plotting of shapes
Dil_growthrates.df <- Dil_growthrates.df %>% arrange(Inoculum, Species, Temp)
# plot all the data together
ggplot(Dil_growthrates.df,
aes(x=Temp, y=mu, colour=Species, group=Species, linetype=Species_in_Expt2)) +
facet_grid(. ~ Inoculum) +
# add a smoothed line to join the different temperatures
#geom_line(aes(x=Temp), stat="smooth", method="loess") +
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=1)) +
geom_point(alpha=0.25, na.rm=TRUE, position=position_dodge(width=1),
shape=Dil_growthrates.df$Shape, size=3) +
geom_linerange(aes(ymin=mu_lwrCI, ymax=mu_uprCI),
alpha=0.25, na.rm=TRUE, position=position_dodge(width=1),
linetype="solid") +
scale_colour_manual(values=ALLsp_palette) +
scale_linetype_manual(values=c("dashed", "solid"))+
labs(x=paste0("Temperature (", "\u00B0", "C)"),
y=expression(paste("Intrinsic Growth Rate (", hr^{-1},")"))
) +
theme(strip.background = element_rect(fill = "white", colour="black"),
strip.text = element_text(color = "black"))## Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
## ℹ Please use the `fun` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
print("Exponential phase culture data below:")## [1] "Exponential phase culture data below:"# plot the exponential data only
muEXP.df <- Dil_growthrates.df %>% filter(Inoculum == "Exponential")
ggplot(muEXP.df,
aes(x=Temp, y=mu, colour=Species, group=Species, linetype=Species_in_Expt2)) +
# do NOT smooth the line used to join species across temperatures
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=1)) +
geom_point(alpha=0.25, na.rm=TRUE, position=position_dodge(width = 1),
shape=muEXP.df$Shape, size=3) +
geom_linerange(aes(ymin=mu_lwrCI, ymax=mu_uprCI),
alpha=0.25, na.rm=TRUE, position=position_dodge(width = 1),
linetype="solid") +
scale_colour_manual(values=ALLsp_palette) +
scale_linetype_manual(values=c("dashed", "solid"))+
labs(x=paste0("Temperature (", "\u00B0", "C)"),
y=expression(paste("Intrinsic Growth Rate (", hr^{-1},")"))
)
# I SHOULD BE USING position_dodge2 TO SHOW ALL THE POINTS,
# but when I do that, the skulls end up in the wrong place
# on the other hand, there are no non-zero small values so maybe we can forget about the shape
ggplot(muEXP.df,
aes(x=Temp, y=mu, colour=Species, group=Species, linetype=Species_in_Expt2)) +
# do NOT smooth the line used to join species across temperatures
stat_summary(alpha=0.5, fun.y=mean, geom="line", position=position_dodge2(width = 1.5)) +
geom_point(alpha=0.25, na.rm=TRUE, position=position_dodge2(width = 1.5)) +
geom_linerange(aes(ymin=mu_lwrCI, ymax=mu_uprCI),
alpha=0.5, na.rm=TRUE, position=position_dodge2(width = 1.5),
linetype="solid") +
scale_colour_manual(values=ALLsp_palette) +
scale_linetype_manual(values=c("dashed", "solid"))+
labs(x=paste0("Temperature (", "\u00B0", "C)"),
y=expression(paste("Intrinsic Growth Rate (", hr^{-1},")"))
)
# plot the stationary data only
print("Stationary phase culture data below:")## [1] "Stationary phase culture data below:"muSTA.df <- Dil_growthrates.df %>% filter(Inoculum == "Stationary")
ggplot(muSTA.df,
aes(x=Temp, y=mu, colour=Species, group=Species, linetype=Species_in_Expt2)) +
# do NOT smooth the line used to join species across temperatures
stat_summary(fun.y=mean, geom="line", position=position_dodge(width=1)) +
geom_point(alpha=0.25, na.rm=TRUE, position=position_dodge(width = 1),
shape=muSTA.df$Shape, size=3) +
geom_linerange(aes(ymin=mu_lwrCI, ymax=mu_uprCI),
alpha=0.25, na.rm=TRUE, position=position_dodge(width = 1),
linetype="solid") +
scale_colour_manual(values=ALLsp_palette) +
scale_linetype_manual(values=c("dashed", "solid"))+
labs(x=paste0("Temperature (", "\u00B0", "C)"),
y=expression(paste("Intrinsic Growth Rate (", hr^{-1},")"))
)
# and plotted again with position_dodge2
ggplot(muSTA.df,
aes(x=Temp, y=mu, colour=Species, group=Species, linetype=Species_in_Expt2)) +
# do NOT smooth the line used to join species across temperatures
stat_summary(alpha=0.5, fun.y=mean, geom="line", position=position_dodge2(width = 1.5)) +
geom_point(alpha=0.25, na.rm=TRUE, position=position_dodge2(width = 1.5)) +
geom_linerange(aes(ymin=mu_lwrCI, ymax=mu_uprCI),
alpha=0.5, na.rm=TRUE, position=position_dodge2(width = 1.5),
linetype="solid") +
scale_colour_manual(values=ALLsp_palette) +
scale_linetype_manual(values=c("dashed", "solid"))+
labs(x=paste0("Temperature (", "\u00B0", "C)"),
y=expression(paste("Intrinsic Growth Rate (", hr^{-1},")"))
)
rm(muSTA.df, muEXP.df)Finally, we can test whether the species’ rank orders remain the same across temperatures. We can do this for all 6 species across all 4 temperatures (using the Kendall’s Coefficient of Concordance). And we can also do it across a pair of 2 temperatures (using Spearman’s Rank Correlation Coefficient).
# get the species rank by finding the median mu of each species
ranks.df <- Dil_growthrates.df %>% group_by(Inoculum, Species, Temp) %>%
summarise(medMu = median(mu))## `summarise()` has grouped output by 'Inoculum', 'Species'. You can override
## using the `.groups` argument.# first try for exponential phase
exp_rank.df <- ranks.df %>% filter(Inoculum == "Exponential") %>%
pivot_wider(names_from = Temp, values_from = medMu)
KendallW(exp_rank.df[,-1*(1:2)], correct = TRUE, test = TRUE)##
## Kendall's coefficient of concordance Wt
##
## data: exp_rank.df[, -1 * (1:2)]
## Kendall chi-squared = 11.385, df = 5, subjects = 6, raters = 4, p-value
## = 0.04427
## alternative hypothesis: Wt is greater 0
## sample estimates:
## Wt
## 0.5692308# we can also test if there's a correlation between 30C and 40C
cor.test(exp_rank.df$`30`, exp_rank.df$`40`, method = "spearman")## Warning in cor.test.default(exp_rank.df$`30`, exp_rank.df$`40`, method =
## "spearman"): Cannot compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: exp_rank.df$`30` and exp_rank.df$`40`
## S = 24.351, p-value = 0.5577
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.3042555# hmm, that's not very meaningful as there is a lot of death. Try 30C & 35C
cor.test(exp_rank.df$`30`, exp_rank.df$`35`, method = "spearman")##
## Spearman's rank correlation rho
##
## data: exp_rank.df$`30` and exp_rank.df$`35`
## S = 8, p-value = 0.1028
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.7714286# also for sanity-check, let's try 25C and 30C
cor.test(exp_rank.df$`30`, exp_rank.df$`25`, method = "spearman")##
## Spearman's rank correlation rho
##
## data: exp_rank.df$`30` and exp_rank.df$`25`
## S = 8, p-value = 0.1028
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.7714286# also try for stationary phase
sta_rank.df <- ranks.df %>% filter(Inoculum == "Stationary") %>%
pivot_wider(names_from = Temp, values_from = medMu)
KendallW(sta_rank.df[,-1*(1:2)], correct = TRUE, test = TRUE)##
## Kendall's coefficient of concordance Wt
##
## data: sta_rank.df[, -1 * (1:2)]
## Kendall chi-squared = 12.615, df = 5, subjects = 6, raters = 4, p-value
## = 0.02726
## alternative hypothesis: Wt is greater 0
## sample estimates:
## Wt
## 0.6307692# test if there's a correlation between 30C and 40C
cor.test(sta_rank.df$`30`, exp_rank.df$`40`, method = "spearman")## Warning in cor.test.default(sta_rank.df$`30`, exp_rank.df$`40`, method =
## "spearman"): Cannot compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: sta_rank.df$`30` and exp_rank.df$`40`
## S = 24.351, p-value = 0.5577
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.3042555# and between 30C and 35C
cor.test(sta_rank.df$`30`, sta_rank.df$`35`, method = "spearman")##
## Spearman's rank correlation rho
##
## data: sta_rank.df$`30` and sta_rank.df$`35`
## S = 8, p-value = 0.1028
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.7714286# finally, for sanity-check: 25C and 30C
cor.test(sta_rank.df$`30`, sta_rank.df$`25`, method = "spearman")##
## Spearman's rank correlation rho
##
## data: sta_rank.df$`30` and sta_rank.df$`25`
## S = 8, p-value = 0.1028
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.7714286################################
# just the subset of 4 species
################################
# now we consider just the subset of 4 species
rank4.df <- ranks.df %>% filter(Species %in% c("P. putida", "P. protegens", "P. grimontii", "P. veronii"))
# first try for exponential phase
exp4_rank.df <- rank4.df %>% filter(Inoculum == "Exponential") %>%
pivot_wider(names_from = Temp, values_from = medMu)
KendallW(exp4_rank.df[,-1*(1:2)], correct = TRUE, test = TRUE)##
## Kendall's coefficient of concordance Wt
##
## data: exp4_rank.df[, -1 * (1:2)]
## Kendall chi-squared = 11, df = 3, subjects = 4, raters = 4, p-value =
## 0.01173
## alternative hypothesis: Wt is greater 0
## sample estimates:
## Wt
## 0.9166667# we can also test if there's a correlation between 30C and 40C
cor.test(exp4_rank.df$`30`, exp4_rank.df$`40`, method = "spearman")## Warning in cor.test.default(exp4_rank.df$`30`, exp4_rank.df$`40`, method =
## "spearman"): Cannot compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: exp4_rank.df$`30` and exp4_rank.df$`40`
## S = 2.254, p-value = 0.2254
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.7745967# and between 30C and 35C
cor.test(exp4_rank.df$`30`, exp4_rank.df$`35`, method = "spearman")##
## Spearman's rank correlation rho
##
## data: exp4_rank.df$`30` and exp4_rank.df$`35`
## S = 0, p-value = 0.08333
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 1# also try for stationary phase
sta4_rank.df <- rank4.df %>% filter(Inoculum == "Stationary") %>%
pivot_wider(names_from = Temp, values_from = medMu)
KendallW(sta4_rank.df[,-1*(1:2)], correct = TRUE, test = TRUE)##
## Kendall's coefficient of concordance Wt
##
## data: sta4_rank.df[, -1 * (1:2)]
## Kendall chi-squared = 11, df = 3, subjects = 4, raters = 4, p-value =
## 0.01173
## alternative hypothesis: Wt is greater 0
## sample estimates:
## Wt
## 0.9166667# test if there's a correlation between 30C and 40C
cor.test(sta4_rank.df$`30`, exp4_rank.df$`40`, method = "spearman")## Warning in cor.test.default(sta4_rank.df$`30`, exp4_rank.df$`40`, method =
## "spearman"): Cannot compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: sta4_rank.df$`30` and exp4_rank.df$`40`
## S = 2.254, p-value = 0.2254
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.7745967# and between 30C and 35C
cor.test(sta4_rank.df$`30`, sta4_rank.df$`35`, method = "spearman")##
## Spearman's rank correlation rho
##
## data: sta4_rank.df$`30` and sta4_rank.df$`35`
## S = 0, p-value = 0.08333
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 1Okay, so there is a pretty good correlation in the species rankings across all the temperatures. For the complete set of all 6 species there is a significant positive correlation in the rankings across all 4 temperatures (both exponential and stationary inocula at p<0.05 by Kendall’s coefficient of concordance Wt). If we consider the sub-set of 4 species used in Expt2, then the p-value becomes slightly smaller but still just p<0.05 so it’s nothing to write home about.
If we consider any particular pair of temperatures (using Spearman’s rank correlation rho), it seems that the only thing happening is that we lose a lot of power when we compare 40C to any of the other temperatures. This is because most species go extinct at 40C and this generates a lot of ties. As a result, the p-value becomes non-significant. This is nothing to write home about either.